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Abstract 


The properties of the observed Higgs boson with mass around 125 GeV can be 
affected in a variety of ways by new physics beyond the Standard Model (SM). 
The wealth of experimental results, targeting the different combinations for the 
production and decay of a Higgs boson, makes it a non-trivial task to assess the 
compatibility of a non-SM-like Higgs boson with all available results. In this 
paper we present Lilith, a new public tool for constraining new physics from 
signal strength measurements performed at the LHC and the Tevatron. Lilith 
is a Python library that can also be used in C and C++/R00T programs. The 
Higgs likelihood is based on experimental results stored in an easily extensible 
XML database, and is evaluated from the user input, given in XML format in terms 
of reduced couplings or signal strengths. The results of Lilith can be used to 
constrain a wide class of new physics scenarios. 
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1 Introduction 


The discovery of a Higgs boson with properties compatible with those of the SM and mass 
around 125 GeV at CERN’s Large Hadron Collider (LHC) [1,2] was a major breakthrough. 
Indeed, the Higgs boson was the last elementary particle predicted by the SM remaining to 
be observed. But, more importantly, the Higgs field has a key role in the SM as it triggers the 
breaking of the electroweak symmetry and gives masses to the elementary particles. Precision 
measurements of the properties of the observed boson are of utmost importance to assess its 
role in the breaking of the electroweak symmetry. They could reveal a more complicated Higgs 
sector, indicating the presence of more elementary scalars or compositeness of the observed 
particle, and could also shed light on a large variety of beyond-the-SM (BSM) particles that 
couple to the Higgs boson. Conversely, precision measurements can be used to rule out new 
physics scenarios affecting the properties of the Higgs boson. 

That the mass of the observed Higgs boson is about 125 GeV is a fortunate coincidence as 
many decay modes of the SM Higgs boson are accessible with a modest integrated luminosity 
at the LHC [3]. Hence, complementary information on the properties of the Higgs boson 
were already obtained from the measurements performed during Run I of the LHC at 7- 
8 TeV center-of-mass energy [4,5]. A large variety of models of new physics (both effective 
and explicit ones) can be constrained from the measurements presented in terms of signal 
strengths. These results were used in a large number of phenomenological studies in the past 
three years (see Refs. [6-31] for a sample of recent studies based on the full data collected at 
Run I). 

However, it is not straightforward to put constraints on new physics from the measured signal 
strengths. Indeed, a large number of analyses have already been performed by the ATLAS 
and CMS collaborations. They usually include several event categories, and present signal 
strength results in different ways. Extracting all necessary information from the figures of 
the various publications is a tedious and lengthy task. Moreover, as the full statistical models 
used by the experimental collaborations are not public, a number of assumptions need to be 
made for constructing a likelihood. The validity of these approximations should be assessed 
from a comparison with the results provided by ATLAS and CMS. 

In order to put constraints on new physics from the LHC Higgs results, many groups have 
been developing private codes. Moreover, recently a public tool, HiggsSignals [32], became 
available. HiggsSignals is a FORTRAN code that uses the signal strengths for individual 
measurements, taking into account the associated efficiencies. In this paper, we present a 
new public tool, Lilith. 1 Lilith is a library written in Python, that can easily be used 
in any Python script as well as in C and C++/R00T codes, and for which we also provide a 
command-line interface. It follows a different approach than HiggsSignals in that it uses as 
a primary input results in which the fundamental production and decay modes are unfolded 
from experimental categories. The experimental results are stored in XML files, making it easy 

1 Lilith is a mythological figure often associated with a female demon. It also stands for “light likelihood 
fit for the Higgs”. 
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to modify and extend. The user input can be given in terms of reduced couplings or signal 
strengths for one or multiple Higgs states, and is also specified in an XML format. 

In Section 2, we present the signal strength framework used to encode deviations from the 
SM at the LHC, as well as the experimental results that we use as input in Lilith. The 
parametrization of new physics effects on the observed Higgs boson, as well as derivation of 
signal strengths, are presented in Section 3. All technical details on how to use Lilith and 
the XML formats that we use are then given in Section 4. Constraints derived from Lilith are 
validated in Section 5, and two concrete examples of its capabilities are given in Section 6 . 
Finally, prospects for Run II of the LHC are discussed in Section 7, and conclusions are given 
in Section 8 . 


2 From experimental results to likelihood functions 

2.1 Signal strength measurements 

Thanks to the excellent operation of the LHC and to the wealth of accessible final states 
for a 125 GeV SM-like Higgs boson, the properties of the observed Higgs boson have been 
measured with unforeseeable precision by the ATLAS and CMS collaborations already during 
Run I of the LHC at 7-8 TeV center-of-mass energy [4,5]. LHC searches are targeting the 
different combinations for the production and decay modes of a Higgs boson. The SM 
Higgs boson has five main production mechanisms at a hadron collider: gluon fusion (ggH), 
vector-boson fusion (VBF), associated production with an electroweak gauge boson (WH and 
ZH, collectively denoted as VH) and associated production with a pair of top quarks (ttH ). 2 
Observation of these production modes constrains the couplings of the Higgs to vector bosons 
(VBF, VH) and to third-generation quarks (ggH, ttH). The main decay modes accessible at 
the LHC are H -> 77 , H -> ZZ* -> M, H WW* 2£2i/, H bb and H rr (with 
t = e,/i). They can provide complementary information on the couplings of the Higgs to 
vector bosons (from the decay into ZZ*, WW*, and 77 ) and to third-generation fermions 
(from the decay into bb, rr, and 77 ). Being loop-induced processes, gg —$■ H and H —» 77 also 
have sensitivity to BSM colored particles and BSM electrically-charged particles, respectively. 

The results of the Higgs searches at the LHC are given in terms of signal strengths, / 1 , which 
scale the number of signal events expected for the SM Higgs, n s . For a given set of selection 
criteria, the expected number of events is therefore fi ■ n s + rib, where rif, is the expected 
number of background events, so that /i = 0 corresponds to the no-Higgs scenario and fi — 1 
to a SM-like Higgs. Equivalently, signal strengths can be expressed as 

cr x A x e m 

M = [a x A x e] SM ’ ^ 

2 Current searches do not constrain the associated production with a pair of bottom quarks, whose SM 
cross section is small and which is plagued with the very large QCD background. 
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where A x e is the product of the acceptance and of the efficiency of the selection criteria. 
Two assumptions can subsequently be made: first, the signal is a sum of processes that exist 
for a 125 GeV SM Higgs boson, i.e. a = J2xy a (X) B(H —» Y) for the various production 
modes X 6 (ggH, VBF, VH, ttH) and decay modes Y <G ( 77 , ZZ*, WW*, bb, rr, ...). 
Second, the acceptance times efficiency is identical to the SM one for all processes, that is 
{A x e)x,Y = [(A x e)x,y] SM for every X and Y. These conditions require in particular that 
no new production mechanism (such as pp —> A —» ZH, where A is a CP-odd Higgs boson) 
exist, and that the structure of the couplings of the Higgs boson to SM particules is as in the 
SM. Under these conditions, signal strengths read 

= T,xA A ^ £ )xMX)B{H^Y) = a(X)B(H^Y) 

^ Ex,y( a x £)x,y<x SM (X) B sm (H -)• Y) ^ X,Y <t sm (X) B sm (H -)• Y) ’ U 

where the eff.Y,y are “reduced efficiencies”, corresponding to the relative contribution of 
each combination for the production and decay of a Higgs boson to the signal. These can be 
estimated from the Axe obtained in a Monte Carlo simulation of individual processes. In the 
case of an inclusive search targeting a given decay mode Y (i.e. VX, (A x e)x,y = (Ax e)y), 
eff'y is equal to the ratio of SM cross sections, crx M /(E x cr| M ). 

The signal strength framework used by the ATLAS and CMS collaborations is based on the 
general form of Eq. (2), hence on the assumption that new physics results only in the scaling 
of SM Higgs processes. This makes it possible to combine the information from various 
Higgs searches and assess the compatibility of given scalings of SM production and/or decay 
processes from a global fit to the Higgs data. This framework is very powerful as it can 
be used to constrain a wide variety of new physics models (some examples can be found in 
Ref. [33]). This is the approach that we will follow in Lilith. However, in order to derive 
constraints on new physics, one first needs to construct a likelihood function from the signal 
strength information given in the experimental publications. In particular, combining the 
results from several Higgs searches is non-trivial and deserves scrutiny. 

2.2 Event categories versus unfolded production and decay modes 

The searches for the Higgs boson performed by the ATLAS and CMS collaborations are 
divided into individual analyses usually focusing on a single decay mode. Within each analysis 
several event categories are then considered. Among other reasons, these are designed to 
optimize the sensitivity to the different production mechanisms of the SM Higgs boson (hence, 
they have different reduced efficiencies ef Tx,y)- In order to put constraints on new physics 
from the results in a given event category, one needs to extract the measurement of the signal 
strength and the relevant e&x,Y information from the experimental publication. For example, 
results of the CMS H —y 77 analysis [34], in terms of signal strengths for all categories, are 
shown on the left panel of Fig. 1. With the addition of the reduced efficiencies effx, 77 , also 
given in Ref. [34], combinations of cr(X) B(H —?• 77 ) can be constrained. 
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Figure 1: Signal strength measurements by the ATLAS and CMS collaborations. On the left 
panel, results of the CMS search H —» 77 [34] category per category. On the right panel, 
2-dimensional ATLAS results in which the fundamental production modes are unfolded from 
experimental categories, in the plane (/i(ggH + ttH, Y), /i(VBF + VH, V)) for Y G ( 77 , ZZ*, 
WW*, tt) [4], 


However, several problems arise when constructing a likelihood. First of all, as can be seen 
on the left panel of Fig. 1, only two pieces of information are given: the best fit to the 
data, that will be denoted as (1 in the following, and the 68 % confidence level (CL) interval 
or la interval. The full likelihood function category per category is never provided by the 
experimental collaborations. Assuming that the measurements are approximately Gaussian, 
it is however possible to reconstruct a simple likelihood, L(/r), from this information. In that 
case, —2 log L(/i) follows a y 2 law. From the boundaries of the 68 % CL interval, left and right 
uncertainties at 68 % CL, A/i~ and A/r + , with respect to the best fit point can be derived. 
The likelihood can then be defined as 


\ ( a #) 2 if 

-2 log L(fi) = l (3) 

[ (&)* if + >A, 

with A= A /j, + in the Gaussian regime. While this is often a valid approximation to the 
likelihood, it should be pointed out that signal strength measurements are not necessarily 
Gaussian, depending in particular on the size of the event sample. 

Barring this limitation, Eq. (3) can be used to constrain new physics. However, it requires 
that at least the 68 % CL interval and the relevant reduced efficiencies effyy are provided 
by the experimental collaboration for every individual category. This is very often, but 
not always, the case. Categories are sometimes defined without giving the corresponding 
signal efficiencies (as in, e.g., the CMS ttH analysis [35]), and/or the result is given for 
a (set of) combined signal strength(s) but not in terms of signal strengths category per 
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category (as in the ATLAS ZZ * and rr analyses [36,37] and in the CMS ttH analysis [35]). 
Such combined \i should in general not be used because they have been obtained under the 
assumption of SM-like production or decay of the Higgs boson. Whenever the effx.y are not 
given in the experimental publications it is in principe possible to obtain estimates from a 
reproduction of the selection criteria applied on signal samples generated by Monte Carlo 
simulation. However, this turns out to be a very difficult or impossible task. Indeed, searches 
for the Higgs boson typically rely on complex search strategies to optimize the sensitivity, 
such as multivariate analysis techniques that are impossible to reproduce in practice with 
the information currently available. Whenever the information on reduced efficiencies is not 
available, we are left to guesswork, with a natural default choice being that e&x is equal to 
the ratio of SM cross sections, d| M /(XA cr v M ), which would correspond to a fully inclusive 
search. 

Constraining new physics from a single LHC Higgs category can already be a non-trivial 
task and come with some uncertainty because the full information is not provided category 
per category. However, more severe complications typically arise when using several cate¬ 
gories/searches at the same time, as is needed for a global fit to the Higgs data. The simplest 
solution is to define the full likelihood as the product of individual likelihoods, 

lw-UUh) ■» -£**(»)w 

i=1 i=l i= 1 \ / 

However, this assumes that all measurements are completely independent. We know that 
this is not the case as the various individual measurements share common systematic uncer¬ 
tainties. They are divided into two categories: the shared experimental uncertainties, coming 
from the presence of the same final state objects and from the estimation of the luminosity, 
and the shared theoretical uncertainties, dominated by the contributions from identical pro¬ 
duction and/or decay modes to the expected Higgs signal in different categories [38]. The 
estimation of the experimental uncertainties in ATLAS should be largely independent from 
the one in CMS, hence these correlations can be treated separately for measurements per¬ 
formed by one collaboration or the other. Conversely, theoretical uncertainties are estimated 
in the same way in ATLAS and CMS and should be correlated between all measurements. 

In the case where all measurements are well within the Gaussian regime, it is possible to take 
these correlations into account in a simple way, defining our likelihood as 

- 2 log L(n) = y 2 (/i) = (/i- (i) T C~\fx - jx ), (5) 

where C is the inverse of the n x n covariance matrix, with Cij = cov[/q, fij] (leading to 
Ca = of). However, the off-diagonal elements of this matrix are not given by the experi¬ 
mental collaborations and are very difficult to estimate from outside the collaboration. This 
remarkably simple and compact expression for the likelihood (a n x n matrix) is only valid 
under the Gaussian approximation; beyond that the expression and the communication of 
the likelihood become more complicated. 
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An alternative way for constraining new physics from the experimental results is to consider 
results in which the fundamental production and decay modes are unfolded from experimental 
categories. These so-called “signal strengths in the theory plane” are defined as 


lAX ,Y) 


o{X) B(H —f Y) 
<7 SM (A')B SM (ff -> Y) ’ 


( 6 ) 


where as before X labels the production mode and Y the decay mode of the Higgs boson. 
These quantities can be estimated from a fit to the results in several event categories; as 
the eff x,y will differ from measurement to measurement, complementary information on the 
various (A", Y ) couples can be obtained and break the degeneracies. The resulting signal 
strengths are directly comparable to the predictions in a given new physics model. They 
have first been used in phenomenological studies in Refs. [39,40]. 

It has become a common practice of the ATLAS and CMS collaborations to present such re¬ 
sults in 2-dimensional likelihood planes for every decay mode. In that case, the five production 
modes of the SM Higgs boson are usually combined to form just two effective A" modes, VBF 
+ VH and ggH + ttH. The likelihood is then shown in the (/i(ggH + ttH, Y), /i(VBF + VH, Y)) 
plane. The ATLAS results in this 2-dimensional plane for Y G ( 77 , ZZ*, WW*, rr), as given 
in Ref. [4], are shown on the right panel of Fig. 1 (the corresponding CMS results can be 
found in Fig. 5 of Ref. [5]). The solid and dashed contours delinate the 68 % and 95% CL 
allowed regions, respectively. As the unfolding of the individual measurements is done by 
the experimental collaborations themselves, all correlations between systematic uncertainties 
(both experimental and theoretical) are taken into account for a given decay mode Y, and 
are encompassed in the correlation between //(ggH + ttH, Y) and //(VBF + VH, Y). (Other 
2-dimensional planes can be relevant, depending on the sensitivity of the searches.) This is a 
very significant improvement over the naive combination of categories of Eq. (4), in which all 
measurements are assumed to be independent. Moreover, in this approach no approximation 
needs to be made because of missing information on the signal efficiencies or signal strengths 
category per category. For these reasons, we use the results in terms of signal strengths in 
the theory plane as the primary experimental input in Lilith. 

A remark is in order regarding the grouping of the five production modes into just two. First 
of all, grouping together VBF, WH and ZH is unproblematic for testing the vast majority 
of the new physics models because custodial symmetry requires that the couplings of the 
Higgs to W and Z bosons scale in the same way. Probing models that violate custodial 
symmetry based on this input and on the inclusive breaking into the individual production 
modes VBF, WH, and ZH, as will be explicited in Eq. (17), may lead to results that deviate 
significantly from the ones using the full likelihood, as will be shown in Section 5.2. The 
combination of the ggH and ttH production modes is more problematic at first sight. While 
gluon fusion is dominated by the top-quark contribution in the SM, this can be modified 
drastically if BSM colored particles are present. However, for all decay modes except H —?• bb 
(where gluon fusion-initiated production of the Higgs is not accessible) the ttH production 
mode is currently constrained with much poorer precision than ggH because of its small cross 
section (being 150 times smaller than ggH at y/s = 8 TeV [3]). Therefore, with the current 



data it is justified to take /i(ggH + ttH, Y) = p(ggH, Y) for all channels except H —> bb , and 
/i(ggH + ttH, bb) = /i(ttH, bb). 3 

Finally, note that all results given in terms of signal strengths are derived assuming the 
current theoretical uncertainties in the SM predictions. Hence, constraining a scenario with 
different (usually larger) uncertainties from a fit to the signal strength measurements is a 
delicate task. This issue will also be discussed, alongside with a possible solution, in Section 7. 

2.3 Statistical procedure 

We use signal strengths for pure production and decay modes as basic ingredients for the 
construction of the Higgs likelihood in Lilith. However, the full likelihood in the /i(X, Y) 
basis is not accessible as only 1- and 2-dimensional (ID and 2D) results are provided by the 
experimental collaborations; therefore some of the correlations are necessarily missing. In 
the currently available ID and 2D results, the full likelihood is provided in some cases in 
addition to contours of constant likelihood. This is extremely helpful since the transmission 
of the result between the collaboration and the reader does not cause any loss of information. 
Two examples from the CMS collaboration are given in Fig. 2. The ID likelihood as a 
function of g(VH, bb) [41] is shown on the left panel . 4 On the right panel, the full likelihood 
in the 2D plane (gt(ggH + ttH, 77 ), /i(VBF + VH, 77 )) [34] is shown as a “temperature plot”. 
Moreover, likelihood grids have been provided by ATLAS in numerical format in the 2D 
plane (/i(ggH + ttH, Y), /i(VBF + VH, Y)) for Y G ( 77 , ZZ\ WW*) [43-45], 

Any result given in terms of signal strengths can be used in Lilith. Whenever available, 
we take into account the full likelihood information. The provision of numerical grids for 
the di-boson final states by the ATLAS collaboration was an important step forward in the 
communication of the likelihood. Unfortunately, they were derived with previous versions of 
the analyses, and the same information has not (yet) been given for the corresponding final 
Run I results [36,46,47]. Moreover, in the CMS H —> 77 result shown on the right panel of 
Fig. 2, the Higgs boson mass has been profiled over instead of being fixed to a given value, 
making the interpretation of the result very difficult. Limitations of current way of presenting 
signal strength results, as well as possible improvements, will be discussed in Section 7. 

If only contours of constant likelihood (the 68 % CL interval in ID, 68 % and 95% CL contours 
in 2D) are present, assumptions about the shape of the likelihood have to be made in order 
to reconstruct it in the full plane. The ID case was already discussed above, and resulted 
in the likelihood of Eq. (3). In the 2D case, a natural choice is to use a bivariate normal 
(Gaussian) distribution. For two (combination of) production and decay processes (A", Y) 

3 Constraints on the ttH production mode for decay modes other than bb are taken into account indepen¬ 
dently in Lilith, see Table 1. 

4 Note that 2D results in the plane (/i(WH, &6), p(ZH, bb)) also exist for this analysis [42]. Both results are 
present in the database of Lilith; by default we use the ID results shown on the left panel of Fig. 2. 
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Figure 2: Signal strength results from the CMS collaboration: ID likelihood 
for VH, H —>• bb (red curve) [41] (left), and temperature plot in the plane 
(MggH + ttH, 77), /i(VBF + VH, 77)) [34] (right). 


and (X', Y'), we obtain the following likelihood: 

-2 log L{n) = ft), 


(7) 


where /1 = 


A X ’> Y '), 


and C 1 = 



is the inverse of the covariance matrix. Under 


the bivariate normal approximation, the 68% and 95% CL contours (which are iso-contours 
of —21ogL) are ellipses. The information on a single contour suffices to reconstruct the 
likelihood in the full plane: the parameters a, b and c, as well as p>{X, Y ) and jl(X\ Y'), can 
be fitted from points sitting on the 68% CL or 95% CL contours as they have known values of 
—21ogL (2.30 and 5.99, respectively). In the following, unless stated otherwise, we choose to 
reconstruct the full likelihood from a fit to the 68% CL contour provided by the experimental 
collaboration. However, having more than one contour of constant likelihood is very useful 
for checking the validity of this approximation. This will be presented in Section 5 for the 
experimental results included in the database of Lilith. Finally, note that generalization of 
the previous equations is trivial should higher-dimensional signal strength measurements be 
published by the experimental collaborations. 


A database of up-to-date experimental results is shipped with Lilith, along with recom¬ 
mended sets of results to use for computing the likelihood (in the form of list hies; all techni¬ 
cal details will be given in Section 4). The default set of results, latest.list, includes the 
latest measurements from the LHC. Its content as of February 2015 is displayed in Table 1. 
All considered 2D results are in the plane (/i(ggH, V), /i(VBF + VH, V)) except for Y = 77 
in ATLAS, where only VBF is considered instead of VBF+VH, and for Y = bb in CMS, 
which is given in the plane (/i(ttH, bb), /i(VH, bb)). The CMS 2D results are taken from the 
combination of Ref. [5], but correspond to the results from Refs. [34,35,53-55]. We also take 
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Collaboration 

Analysis 

Type 

Reference 


H > 77 

2D contour 

[46] 


H ->■ ZZ* 

2D contour 

[36] 


H ->• WW* 

2D contour 

[47] 

ATLAS 

H —>• rr 

2D contour 

[37] 




VH, H -> bb 

2D contour 

[48] 


ZH, H ->• invisible 

full ID 

[49] 


ttH, H -> bb 

ID interval 

[50] 


ttH, H — » 77 

ID interval 

[46] 


H ->• 77,ZZ*,fFfF*,^,TT 

2D contours 

[5] 

CMS 

ttH, H — > 77, rr 

ID interval 

[5] 


ttH, H — )• leptons 

ID interval 

[35] 


ZH + VBF, H — * invisible 

full ID 

[51] 

CDF & DO 

VH, H ->• 66 

ID interval 

[52] 


Table 1: Recommended set of experimental results, for the database of Lilith version 15.02. 
This set corresponds to the hie latest.list, and is used by default when running Lilith. 








into account all available searches on production in association with a top-quark pair, as well 
as searches for invisible decays of the Higgs boson from both ATLAS and CMS. Note the 
presence of the CDF and DO combined result for VH, H —>• bb [52]; only in this channel the 
precision of the Tevatron result is comparable with the one of the LHC at Run 1. 

All considered experimental results are given at a fixed Higgs mass (that can be read from the 
database, see Section 4.5) in the [125,125.6] GeV range. Variations of the experimental results 
within this narrow interval are expected to be small, hence limiting the inconsistencies when 
combining the results. However, it would be desirable to take into account the variation of 
the results with mass, as we will argue in Section 7. The final Higgs likelihood is the product 
of the individual (1- or 2-dimensional) likelihoods. Validation of the Higgs likelihood used in 
Lilith against official LHC results will be presented in Section 5. 


3 Parametrization of new physics 

In order to assess the compatibility of a new physics hypothesis with the LHC measurements 
presented in the previous section, one needs to compute the expected signal strengths /i(X, Y ) 
(see Eq. (6)) for the relevant production mechanisms X and decay modes Y. This can 
be achieved in a direct way from ct(X), cr SM (X), B(H —> Y), and B SM (H —> Y), but is 
often found to be impracticable. Indeed, in order to have well-defined signal strengths (for 
which /i = 1 corresponds to the SM prediction) one should take the same prescription for 
computing cross sections and branching fractions in the SM and in the considered new physics 
scenario [3]. Concretely, one needs to consider the same order in perturbation theory, the 
same set of parton density functions, etc. 

In most new physics scenarios only leading order (LO) computations are available. Thus, 
all available next-to-LO (NLO) corrections to the SM predictions should be ignored. While 
this leads to properly defined signal strengths, o'NLo(^)/ cr NLo(^) typically differ from 
a LO (X)/(r^(X) (and similarly for branching ratios) as soon as one deviates from the SM 
prediction. This is because the relative contributions of SM particles to the process will be 
affected by the NLO corrections. For instance, higher-order corrections to the gluon fusion 
process will change the relative contribution of the top and bottom quark loops. Therefore, 
considering LO or NLO cross sections will yield different /i(ggH, Y) if new physics affects the 
couplings of the Higgs to top and bottom quarks in a different way. 

These two problems come from the parametrization of new physics effects from cross sections 
and branching ratios. As we will see, they can be alleviated if new physics is parametrized 
instead using reduced couplings. 
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3.1 Scaling factors and reduced couplings 

The general signal strength expression given in Eq. (2) can be rewritten as 


a(X) B(H -» Y) 

** ^ a SM (X) B sm (H Y) 


c 2 x a SM (x) w ^ r| M 
<r SM (x) x r?, M x xvcfrF 


1 

Ey Cl B sm (H Y) 


E effx.yC'iC' 2 , 

X,Y 


( 8 ) 


where Tf ^ 1 is the total decay width of the SM Higgs boson, and the cross section (partial 
width) for each process X (Y) is scaled with a factor C\ (Cy) compared to the SM expec¬ 
tation. 5 The term YIy Cy B SM (H —> Y) accounts for the scaling of the total width of the 
Higgs boson. (We assume that the narrow-width approximation also holds in the new physics 
scenarios.) Furthermore, we introduce reduced couplings through the following Lagrangian, 


£ = 9 
+ 9 


CwrnwWW,, + C z 


m z 


cos 6 


-Z»Z n 


w 


H 


n m t ,7 n rn b ,7 n m c _ m T _ 
—C t - - tt — 6 ft-- bb — C c - - cc — C r --rr 

2m\y 2m\y 2m\y 2m\y 


H. 


(9) 


where Cw,z and C t} b,c,r are bosonic and fermionic reduced couplings, respectively. Light 
fermions are not taken into account as their phenomenological impact on the SM Higgs sector 
is negligeable. In the limit where all reduced couplings go to 1, the SM case is recovered. 
At leading order in perturbation theory, the scaling factors Cx and Cy from Eq. ( 8 ) can be 
directly identified with the reduced couplings Ci from Eq. (9) for processes involving just one 
coupling to the Higgs boson. We obtain 


r^2 _ si 2 s~i2 _ s~y2 s~i2 _ f~~i2 _ _ r^i2 s~^2 _ r^2 

u wh — ? °zh — ? °ttH — > °// — °/ > — > 


( 10 ) 


where / = b, c, r and V — W, Z. 

For the remaining main processes (ggH and VBF production, decay into gg, 77 and Z 7 ), 
there is no direct identification unless the Higgs couplings to all involved SM particles scale 
in the same way. In the general case, the Cx and Cy for these processes will be given by a 
combination of reduced couplings C t , weighted according to the contribution of the particle 
i to the process. For the production mechanisms, we have 


E C,C>if(ggH) 

=t,b,c 

, E ^ M (ggH) 


E CiCj erf M (VBF) 

r 2 _ ij=W,Z ___ 

VBF E <nf(VBF) ’ 

i,j=W,Z 


( 11 ) 


where the <r, s , M are the different contributions to the cross section in the SM. For i = 7 
it corresponds to the cross section from the particle i alone, while j ^ j comes from the 

5 The scaling factors Ci are often seen elsewhere as Ki. 
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interference between the particles i and j. (We only consider either the term or in 
the sum, not both, to avoid double counting the interference terms.) Similarly, the reduced 
couplings for the gg, 77 , and Z 7 loop-induced decay modes are computed as 


£ CtCjTffiHgg) 

— i <3~ t ’ b ’ c _ s~i2 

^ r->• gg) ’ 

i,j=t,b,c 


E CiCj If (II -> 77 . z 7 ) 

i,j=W,t,b,c,T 

e 

i,j=W,t,b,c,T 


( 12 ) 


where the are the SM partial widths of the process under consideration. In all cases, all 
relevant SM contributions have been taken into account. Note that the relative sign of the 
reduced couplings will affect the interference terms, as they are proportional to CiCj. 


At LO, the various and can be obtained from tree-level amplitudes (for VBF) or from 
the 1-loop amplitudes (for gg —?• H and H —» gg, 77 , Z 7). 6 It would however be desirable to 
take into account the NLO corrections to the Higgs cross sections and partial widths as they 
modify the relations Cx,y(Ci). This can be achieved in a simple way as long as higher-order 
corrections only rescale the crf) M and T® M that are already existing in Eqs. (1 1)—(12) , and do 
not induce new couplings to the Higgs boson. This is the case for the QCD corrections, but 
not for the electroweak corrections. Thus, as will be explained in Section 4.6.2, the QCD 
corrections for all Eve processes of Eqs. (11)—(12) will be included in Lilith. 


One last remark is in order. The signal strength framework requires that the signal in all 
searches be a sum of processes that exist for the SM Higgs boson. However, new production 
or decay modes may exist without spoiling the signal strength interpretation as long as they 
do not yield sizable contribution in the current Higgs searches. Two particularly interesting 
cases are Higgs boson decays into undetected particles, or into invisible particles. In the first 
case, this new decay is simply missed by current searches (as would, e.g., be the case for the 
decay of the Higgs into light quarks and gluons), while in the second case this new decay 
mode gives rise to missing energy in the detector. As was shown in Section 2.3, invisible 
decays of the Higgs boson are constrained by current searches which are taken into account 
in Lilith. In both cases of undetected and invisible decays, the width of the Higgs boson 
becomes larger and modifies the signal strength predictions of Eq. ( 8 ) as 

g{Cxi Cy) t (1 I3; nv i s ible ^undetected )fJ>( L Cxi Cy) . (13) 

In Lilith, arbitrary invisible and/or undetected decays can be specified, as will be presented 
in Section 4.6. 


3.2 CP-violating admixtures 

We also consider the case where the observed Higgs boson is a mixture of CP-even and 
CP-odd states [56,57]. The Higgs coupling to vector bosons has the form 

VVH : C v g -^~ g ,lu , (14) 

mw 

6 At LO, one obtains the same scaling for gluon fusion and for the decay into two gluons, C gg h = C gg . 
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where as above Cy measures the departure from the SM: Cy = 1 for a pure scalar H + 
(CP-even) state with SM-like couplings and Cy = 0 for a pure pseudoscalar H~ (CP-odd) 
state. 

In the fermion sector, we find the general vector and axial-vector structure of the Higgs 
coupling to fermions. Concretely, we have 

Hff: + (15) 

Z7Tl\y 


where in the SM one has Re(Cy) = 1 and Im (Cf) = 0, while a purely CP-odd Higgs would 
have Re(C/) = 0 and Im(C/) = 1. Since rnj <C m 2 H for / = 6 , c, r, the partial decay widths 
scale as T(H —> ff) oc Re(C /) 2 + Im(C /) 2 = |C /| 2 to a very good approximation [58]. This 
is what is implement in Lilith. Effects of CP mixing will mainly show up at loop level, in 
particular in the gg —>• H and H —y 77 rates. A test of the CP properties of the observed 
Higgs from a global fit to the signal strengths was presented in [9,59]. Following Ref. [9], at 
leading order the Higgs rates normalized to the SM expectations can be written as 
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(16) 


with A/ [rri\y\ ~ —8.32 for m h = 125 GeV. For convenience, the contribution from the other 
quarks has been omitted in the above equations but is taken into account in Lilith. 

In the case of ttH production, the approximation that we used above for the other fermions 
does not hold since m t > mu- Instead, the cross section scales as cr(ttH +// ~) oc Re(C /) 2 + 
Im(C'/) 2 cr SM (ttH _ )/cr SM (ttH + ). Following Ref. [60], a factor cr SM (ttH~)/cr SM (ttH + ) ss 1/3 
is considered in Lilith. However, a significant coupling of the CP-odd component of the 
Higgs boson to top quarks may modify the acceptance times efficiency compared to the SM 
value in searches for the Higgs boson in association with a pair of top quarks [60,61], i.e., 
(.A x e)ttH,y 7^ [(A x e)ttH,y] SM - As this cannot be taken into account in Lilith, such cases 
should be interpreted with care. Moreover, only after the end of Run II will the LHC have 
enough sensitivity to probe CP violating effects in the H —>• rr decays [62] , and the product 
(A x £)x,tt can thus be approximated by the SM one for now. Details on how to specify real 
and imaginary parts for the couplings are given in Section 4.6.2. More precise measurements 
at Run II of the LHC will ultimately call for an implementation of CP admixture that includes 
NLO effects in Lilith. 
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Running Lilith 


4.1 Getting started 

Lilith is a library written in Python for constraining model of new physics against the 
LHC results. The code is distributed under the GNU General Public License v3.0. The 
latest version of Lilith and of the database of experimental results (as of February 2015, 
Lilith 1.1 and database version 15.02) as well as all necessary information can be found 
at 


http://lpsc.in2p3.fr/projects-th/lilith 

The archive of Lilith can be unpacked in any directory. It contains a root directory called 
Lilith-1.1/ where the following directories can be found: 

• lilith/: the Python package itself. The Lilith application programming interface 
(API) will be presented in Section 4.2. It also contains the Python/C API that will be 
presented in Section 4.4. 

• data/: contains the database of experimental results in XML format, as well as * . list 
text hies for the recommended sets of results. Details are given in Section 4.5. 

• userinput/: where parametrizations of new physics models, in the XML format de¬ 
scribed in Section 4.6, can be stored. Some basic user input hies that include extensive 
comments are provided with the Lilith distribution. 

• examples/: concrete examples on how to use Lilith for constraining new physics. 
Two of them will be presented in detail in Section 6; an example for using Lilith in C 
and C++/R00T programs will be discussed in Section 4.4. 

• results/: empty folder where results from Lilith can be stored. 

The folder Lilith-1.1/ moreover contains run_lilith. py, the command-line interface (CLI) 
of Lilith that will be presented in Section 4.3, as well as general information, information 
on the license, and a changelog in the hies README, COPYING, and changelog, respectively. 

Lilith requires Python 2.6 [63] or more recent, but not the 3.X series. The standard 
Python scientific libraries, SciPy and NumPy [64], should furthermore be installed. We require 
SciPy 0.9.0 or more recent, and NumPy 1.6.1 or more recent. Python, SciPy and NumPy are 
available for the major platforms, including GNU/Linux, Mac OS X, and Microsoft Windows. 

The easiest way to check if all dependencies of Lilith are correctly installed is to try to 
compute the likelihood from an example hie. This can be achieved by typing to the shell 
(with current directory Lilith-1.1/) the command 
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python run_lilith.py userinput/example_couplings.xml 

Everything is correctly installed if basic information as well as the value of the likelihood is 
printed on the screen. Note that the version number of Python can be obtained by typing 
the command python —version to the shell, while the presence of SciPy and NumPy and 
their version numbers can be checked by typing in an interactive session of Python (started 
by typing python to the shell) the following commands: 

import scipy 

print scipy._version_ 

import numpy 

print numpy._version_ 

Note that every version of Lilith is shipped with the latest version for the database of 
experimental results at the time of release. However, as new experimental results usually do 
not require any modification to Lilith, we do not release a new version of the code every time 
new experimental results come out. Instead, we provide separately an update of the database 
of experimental results. Each release of the database has version number YY. MM (e.g., 15.02), 
where YY and MM correspond to the year and to the month, respectively. If two or more updates 
of the experimental database are provided the same month, from the second release onwards 
versions will be numbered YY.MM.n (with n starting from 1). The version number of the 
database can be found in data/version (also accessible via the readdbversion() method 
of the API, see next section, and printed on the screen when using the CLI). When using 
the recommended sets of experimental results of Lilith in a publication, the version number 
of the database must also be cited in addition to the experimental publications from which 
results were used. 


4.2 The Lilith API 

Lilith provides an API from which all tasks (reading the user and the experimental input, 
compute the likelihood, print the results in a hie, etc.) can be performed, using the methods 
described below. This is the recommended way of using Lilith. In order to be used in any 
Python code (or in an interactive session of Python), the package of Lilith, called lilith, 
first needs to be imported. However, Python needs to know the location of the lilith 
package. This can be achieved in at least three ways: 

1. create the Python script importing lilith in the directory Lilith-1.1/ (or, in an 
interactive session, having Lilith-1.1/ as the current directory). 

2. adding the path to Lilith-1.1/ to the environment variable PYTHONPATH. This can be 
done with the command 

i I export PYTHONPATH=/path/to/Lilith-l.l:$PYTHONPATH 
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in bash shell or 

setenv PYTHONPATH /path/to/Lilith-1.1:$PYTHONPATH 


in esh/tesh shell. In order to permanently have the path to Lilith in PYTHONPATH 
(not only for the current session) this command should be added in a .bashre (for 
bashf) or .eshre / .teshre (for esh/tesh shell) hie located in the home directory of 
the user. 

3. adding the path to Lilith-1. 1/ to the variable sys .path by starting the script with 

1 import sys 

2 sys.path.append( 1 /path/to/Lilith-1.1') 

before importing lilith. Note that the path can also be relative. 

The Lilith library can then be imported to the current script by typing import lilith 
or from lilith import *. In the first case, all classes, methods and attributes defined in 
the code will be in the namespace lilith, in the second case they will be in the global 
namespace. We now present all methods and attributes of the API of Lilith. 

class lilith.Lilith (verbose—False, timer=False) 

Instanciate the Lilith class. The following public attributes are initialized: 

verbose 

if True, information will be printed on the screen 

timer 

if True , each operation will be timed and the results will be printed on the screen 

expjnu 

list of experimental results read from the database 

exp_ndf 

number of measurements (n-dimensional results count for n measurements) 

dbversion 

version number for the database of experimental results 

couplings 

list of reduced couplings for each Higgs particle contributing to the signal as read 
from the user input 

user mu 

list of signal strengths for each Higgs particle contributing to the signal as read 
(or derived from) the user input 

user_mu_tot 

signal strengths for the sum of the Higgs particles present in user_mu 
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results 

list of results after computation of the likelihood for each individual measurement 

1 

value of —2 log L 
exception lilith.LilithError 

Base exception of Lilith; all other exceptions inherit from it. For the definition of all 
the exceptions, see lilith/errors .py. 

Lilith.readuserinput (userinput) 

Read the string in XML format given as argument and fill the attribute couplings (if 
the user input is given in terms of reduced couplings) or user_mu and user_mu_tot (if 
the user input is given in terms of signal strengths). User input formats are presented 
in Section 4.6. 

Lilith. readuserinputfile(/f/epat/i) 

Read the user input located at filepath and call readuserinput(). 

Lilith. computecouplings() 

Compute from user_couplings the following reduced couplings if not already present 
in user .couplings: CVbf, CggH, C gg , C 77 , and Czy 

Lilith. computemufromreducedcouplingsQ 

Compute the signal strengths (stored in user_mu and user_mu_tot) from the reduced 
couplings in user .couplings. 

Lilith. compute_user_mu_tot() 

Add up the signals from all Higgs bosons contributing to the signal in user_mu; store 
the result in user mu tot. 

Lilith, readexpinput ( filepath=default_expJist ) 

Read the experimental input specified in a list hie and store the results in exp_mu 
and exp ndf. By default, the list hie is data/latest. list. The formats of the 
experimental results are presented in Section 4.5. 

Lilith. readdbversion() 

Read the version of the database of experimental results from the hie data/version, 
and store the information in dbversion. 

Lilith. compute exp ndf() 

Compute the number of measurements from exp_mu and store the information in 

exp_ndf. 

Lilith. computelikelihood(?zsennjmf=iVone, exp-filepath=None, userfilepath=None ) 

Evaluate the likelihood function from signal strengths derived from the user input 
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(user mu tot) and the experimental results (exp mu) and store the results in the at¬ 
tribute results. If the arguments userinput and userfilepath are not specified, user mu tot 
will be assumed to have been filled already. Else, all information will be read from the 
XML input given in userinput , or from the hie located at userfilepath, and user_mu_tot 
will be computed. If the exp-filepath argument is not specified, the experimental results 
will be read from the default list hie unless exp mu is already filled. Else, experimental 
results from exp-filepath will be read before computing the likelihood. 

Lilith, writecouplings [filepath) 

Write reduced couplings from the attribute couplings in a hie located at filepath in 
the XML format specihed in Section 4.6.2. 

Lilith.writesignalstrengths(/i/epat/i, tot=False ) 

Write signal strengths from the attribute userunu (if tot=False) or user mu tot (if 
tot=True ) at filepath in the XML format specihed in Section 4.6.1. 

Lilith .writeresults (filepath, slha=False) 

Write the content of the attribute results at the location filepath in the XML format 
(if slha—False) or the SUSY Les Houches Accord (SLHA)-like format [65] specihed in 
Section 4.7. 

Note also that the version of Lilith is stored in the lilith. __version__ variable. 

A minimal example of use of the API is as follows: 

1 from lilith import * 

2 leal = LilithO 

3 leal.readexpinput() 

4 leal.readuserinputfile( 1 userinput/example_mu.xml 1 ) 

5 leal. computelikelihoodO 

e print 1 - 21 og(likelihood) =', lcal.l 

The hrst two lines import the Lilith library into the global namespace and initialize the 
computations. They are equivalent to 

1 import lilith 

2 leal = lilith. LilithO 

The three following lines successively read the experimental input, read the user input from 
the hie userinput/examplennu.xml, and compute the likelihood. Alternatively, they could 
be replaced with a single line, 

i | leal.computelikelihood(userfilepath='userinput/example_mu.xml') 

Finally, the value of —21ogL is printed on the screen on the last line. 

In the example above, any error (corresponding to an exception in Python) will interrupt 
the code. This may not be the desired behavior. In particular, if several user inputs are 
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successively given to Lilith (as in the case of a scan of a parameter space), it may be 
preferable to store the error and move on to the next user input instead of stopping the 
execution of the code. In Python, the handling of errors can be achieved with try . . 
except blocks. We provide below a simple example. 

1 try: 

2 leal.readuserinputfile('userinput/example_mu.xml') 

3 leal. computelikelihoodO 

4 print 1 - 21 og(likelihood) =', lcal.l 

5 except LilithError as e: 

e print 'the following error occured: 1 , e 

Here, any error raised by Lilith (of type lilith. LilithError or derived from it) will be catched, 
in which case the error message will be printed on the screen and the script will continue 
normally. It is of course also possible to store the error message in a hie, or simply replace the 
last line with the pass statement in order to ignore all errors raised by Lilith and continue 
the execution of the script. For the definition of all exceptions used in Lilith (that all derive 
from lilith. LilithError), see lilith/errors. py. This makes it possible to treat each type 
of error in a different way. 

4.3 Command-line interface 

A command-line interface or CLI is also shipped with Lilith for a more basic usage of the 
tool. It corresponds to the hie run_lilith.py located in the directory Lilith-1.1/. The 
CLI can be called by typing to the shell (with current directory Lilith-1.1/) the command 

i | python run_lilith.py user_input_file (experimental_input_file) (options) 

where arguments in parentheses are optional. 

The hrst argument, user_input_f ile, is the path to the user input hie in the XML format 
described in Section 4.6. New physics can be parametrized in terms of reduced couplings 
(see Section 3) or directly in terms of signal strengths. Examples are shipped with Lilith in 
the directory userinput/. The second argument, experimental_input_f ile, is the path to 
the list of experimental results to be used for the construction of the likelihood. If not given, 
the latest LHC results will be used (data/latest. list; its content is given in Table 1). It 
is the recommended list of experimental results to be used for performing a global ht. All 
details about the experimental input will be given in Section 4.5. 

If no option is given, basic information as well as the value of the likelihood and the number 
of measurements is printed on the screen. A number of options are provided to control the 
information printed on the screen and to print the results of Lilith in output hies. They 
are listed in Table 2. 
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Option 

Meaning 

—help, -h 

Display basic usage of runJ.ilith.py 
and list of options 

—verbose, -v 

Display details about the computation 

—timer, -t 

Time each operation and display results on the screen 

—silent, -s 

Silent mode 

—couplings=output, 

-c output 

Obtain and print the attribute couplings 

in the hie output in XML format 

—mu=output, 

-m output 

Obtain and print the attribute user_mu 
in the hie output in XML format 

—results=output, 

-r output 

Obtain and print the attribute results 
in the hie output in XML or SLHA-like format 


Table 2: Options available when r unning run_lilith.py. 


The option —couplings / -c only works in the reduced couplings mode. In addition to 
the reduced couplings already present in the input hie, it prints scaling factors computed 
from the input (he. C' 77 , C'^ 7 , C ggi C' gg H, and CVbf, see Section 3.1). The option —mu / -m 
prints the complete list of signal strengths in a hie in XML format. More specifically, all signal 
strengths fi(X,Y) with X e (ggH, VBF, WH, ZH, ttH) and Y e ( 77 , ZZ *, WW*, bb, rr, 
cc, Z 7 , gg, invisible) are printed. Finally, the option —results / -r prints the value of 
—21ogL and the number of measurements in a hie. If the extension of the filename is . slha 
(case-insensitive), a hie in SLHA-like format is created. Otherwise, an XML hie is created, 
with extra information on the individual contributions to —2 log L from all the experimental 
results used in the calculation. More details on the structure and content of the output hies 
are given in Section 4.7. 

Initialization of Lilith and reading of the database of experimental results is done at each 
execution of run_lilith.py. When computing the Higgs constraints in the context of the 
scan of a model, this only needs to be done once. Therefore, successive calls to the CLI will 
be much slower than direct calls to the methods of the API (see previous section), even more 
so as internal information is stored to optimize successive computations of the results when 
using the API. Whenever performance can be an issue, we highly recommend using the API 
instead of the CLI. 
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4.4 Interface to C and C++/R00T 


Above, we have presented how to use Lilith from the API written in Python and from the 
CLI. In order to use Lilith within a C, C++, or ROOT code, a first possibility is to call the CLI. 
However, this will suffer from the performance issues explained at the end of the previous 
section. Fortunately, Python provides a Python/C API [66] from which each method of the 
API of Lilith can be called, and each attribute of Lilith can be manipulated, in both C 
and C++. However, direct use of the functions of the Python/C API can be quite tedious. 
For that reason, we have also written a C API for Lilith, consisting in a series of functions 
using the Python/C API for performing all usual tasks, following closely the methods of the 
API presented in Section 4.2. 

The C API for Lilith is contained in the directory lilith/c-api. Its use requires the 
libraries and header files needed for Python development. They can be installed from most 
package managers under the name python-dev or python-devel, depending on the platform. 
More information on how to install these libraries and header, platform by platform, can be 
found at [67]. 

We now present all the functions of the C API: 
initialize_lilith( char* experimental-input) 

Import lilith, instanciate the class Lilith, and read the experimental input file located 
at experimental-input. The function returns the instance object. If experimental-input 
is an empty string(""), the default experimental input file data/latest.list is used. 

hhthj'eaduserinput(PyOF/eci* lilithcalc, char* XMLinputstring) 

Read the user input XML string XMLinputstring and store the information in the object 
lilithcalc. 

lilith_readuserinput_fromhle(PyO&ject* lilithcalc , char* XMLinputpath) 

Read the user input XML file located at XMLinputpath and store the information in the 
object lilithcalc. 

hlith_computelikelihood(Py(9&ject ,ir lilithcalc) 

Evaluate and return the value of —21ogL from the object lilithcalc. 

lilith_exp_ndf (PyObject* lilithcalc ) 

Evaluate and return the number of measurements from the object lilithcalc. 

hlith_likelihood_output(Py(9&ject* lilithcalc , char* outputfilepath, int slha) 

Write the content of the attribute results at outputfilepath in the XML format (if slha= 0) 
or the SLHA-like format (otherwise) specified in Section 4.7. 

hhth mu output (PyO&ject* lilithcalc , char* outputfilepath , int tot) 

Write signal strengths from the attribute user_mu (of tot=0) or user_mu_tot (other¬ 
wise) at outputfilepath in the XML format specified in Section 4.6.1. 
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lilith couplings output(JA/O&jeci* lilithcalc, char* outputfilepath) 

Write reduced couplings from the attribute couplings in a file located at outputfilepath 
in the XML format specified in Section 4.6.2. 

Initialization of Lilith and the reading of the experimental input hie can be done just once 
with the function initialize lilith(). Evaluation of the likelihood can then be performed 
separately. An example of use of the C interface of Lilith is shipped with the code. It is 
available at examples/c/lilith_compute. c. We now present it step by step. 

#include <Python.h> 

#include "lilith.h" 

int main(int arge, char* argv[]) 

Py_Initialize(); 

The Python.h and lilith.h are the Python/C API and Lilith header hies, respectively. 
Those are linked during the compilation by the Makefile located in the same directory. To 
start the Python/C API, the function Py_Initialize() is mandatory. 

char experimental_input [] = /data/latest.list"; 

char output_couplings[] = "lilith_couplings_output.xml"; 
char output_XML[] = "lilith_likelihood_output.xml"; 
char output_SLHA[] = "lilith_likelihood_output.slha"; 
char output_mu[] = "lilith_mu_output.xml"; 

Line 6 is the path to the the experimental input hie. Note that in this particular example, 
it is equivalent to char experimental.input[] = Various output hie paths are then 
defined. 

PyObject* lilithcalc = initialize_lilith(experimental_input); 

Line 12 is the initialization of a Lilith object lilithcalc from the experimental input hie 
experiment .input. 

char XMLinputstring[6000] = 

[...] // Construction of the user XML input string 

lilith_readuserinput(lilithcalc, XMLinputstring); 

On lines 13 and 14, the user XML string input is constructed, see examples/c/lilith_compute. c 
for more details. The function lilith_readuserinput() is then called to read the user input. 

float my_likelihood; 

my_likelihood = lilith_computelikelihood(lilithcalc); 
printf ("-2*log(L) = °/„lf \n" , my_likelihood) ; 
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int exp_ndf; 

exp_ndf = lilith_exp_ndf(lilithcalc); 
printf ("exp_ndf = 7„i\n", exp_ndf) ; 

lilith_likelihood_output(lilithcalc, output_XML, 0); 
lilith_likelihood_output(lilithcalc, output_SLHA, 1); 
lilith_mu_output(lilithcalc, output_mu, 0); 
lilith_couplings_output(lilithcalc, output_couplings); 

Py_Finalize(); 
return 0; 

} 

Once the experimental and user input files have been read, computations can be performed. 
First, —21ogL and the number of measurements are printed on the screen. Various output 
hies are also created. Finally, Py.Finalize() can be used to free all memory allocated by 
the Python interpreter. 

This example, lilith.compute. c, can be compiled and executed by typing to the shell 

1 make 

2 ./lilith_compute 

while the executable and the intermediate object hies can be removed by typing to the shell 
make clean. In order to use Lilith in C++ or ROOT codes, only mininal changes need to be 
made compared to the case of C codes. C++ users may want to change the compiler (CC in the 
Makefile) from gcc to g++. ROOT users should furthermore link the headers and libraries of 
ROOT by adding the following two lines to the Makefile: 

1 CFLAGS += $(shell root-config —cflags) 

2 LFLAGS += $(shell root-config —glibs) 
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4.5 Experimental input 

We have seen that the evaluation of the likelihood in Lilith requires the input of a list 
of experimental results to be considered. It corresponds to a simple text hie with a .list 
extension listing the paths to experimental result hies in XML format (each containing a 
single ID or 2D signal strength result). Lilith is shipped with the latest LHC Higgs results 
(plus a Tevatron result), see Table 1 in Section 2.3, in the form of XML hies present in 
subdirectories of data/. Moreover, several lists of experimental results are provided in data/, 
with latest .list being the default list hie. This is the one recommended for a global ht to 
the LHC+Tevatron Higgs data. 
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The user can also create his/her own list hie in the data/ directory. For instance, in order to 
put constraints on new physics using only the latest di-boson results from ATLAS [36,46,47], 
one can create a hie list that contains 

# ATLAS di-boson analyses 

ATLAS/Runl/HIGG-2013-08_ggH-VBF_gammagamma_n68.xml 
ATLAS/Runl/HIGG-2013-21_ggH-VVH_ZZ_n68.xml 
ATLAS/Runl/HIGG-2013-13_ggH-VVH_WW_n68.xml 

The first line, starting with a #, is a comment and is not read by Lilith. The three following 
lines indicate the paths to the XML hies to be considered. As can be read in the paths, these 
are published results from the ATLAS collaboration based on Run I data. The conventional, 
though not mandatory, naming scheme for the hies is as follows. The identifier of the analysis 
(HIGG-2013-XX in this case) comes hrst, followed by the 2D plane in which results are given; 
finally, n68 indicates that the likelihood has been reconstructed from the contour at 68% CL 
under the Gaussian approximation. Note that no consistency check is done by Lilith. 
When creating a new list hie, the user should make sure that there is no overlapping between 
experimental results (e.g., that the results of two versions of the same analysis, based on 
overlapping event sets, are not used at the same time). 

4.5.1 XML format 

Every single experimental result (ID or 2D) is stored in a different XML hie. In this way, 
modifying and updating the database is an easy process. We now present the format of the 
experimental results hies in Lilith. 

The root tag of each hie is <expmu>. It has two mandatory attributes, dim and type, that 
specify the type of signal strength result (ID interval, full ID, 2D contour, or full 2D, see 
Section 2.3). Possible values for the attributes are given in Table 3. In addition, the <expmu> 
tag has two optional attributes: prod and decay. They can be given a value listed in Table 3 
if the analysis under consideration is only sensitive to one production mode (e.g., ttH) or 
to one decay mode (e.g., 77 ) of the Higgs boson. In the general case, the prod and decay 
attributes can be skipped. Indeed, all relevant efficiencies ef fx,Y (see Section 2.1) can be 
specihed in <eff> tags. 

Taking for instance the case of a ID measurement, one can specify 

<eff prod="ggH">0.5</eff> 

<eff prod="VBF">0.5</eff> 

if decay is specihed in <expmu>. Similarly, one can specify 

<eff decay="WW">0.5</eff> 

<eff decay="tautau">0.5</eff> 
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Attribute 

ID interval 

full ID 

2D contour 

full 2D 

dim 

II ^ II 

II II 

"2" 

"2" 

type 

"n" 

II II 

" n " 

II II 

(prod) 

"ggH", "ttH", "VBF 

", "WH", "ZH", "VH", "VVH" 

(decay) 

"gammagamma", "ZZ", "WW", "Zgamma", 

"tautau", "bb", "cc", "invisible" 


Table 3: Allowed values for the attributes of the <expmu> tag. in experimental hies in XML 
format. The four different formats of experimental data are defined by the mandatory dim 
and type attributes. 


if prod is specified in <expmu>. If none of them is present, one could specify efficiencies in 
the following way 


<eff prod="ggH" 
<eff prod="VBF" 
<eff prod="ggH" 
<eff prod="VBF" 


decay="WW">0.25</eff> 
decay="WW">0.25</eff> 
decay="tautau">0.25</eff> 
decay="tautau">0.25</eff> 


where it is required that the sum of all efficiencies is 1. In the case of 2D signal strengths (i.e. 
if dim="2" in <expmu>) the efficiencies should be given for both dimensions (and separately 
add up to 1). In this case, the attribute axis should be provided in every <eff> tag, with 
possible values being "x" or "y" for the first and second dimension of the results. 

Before turning to the syntax case by case, we comment on the optional information that can 
be provided in the experimental hies. The following tags can be given: 


• <experiment>, referring to the experiment that has produced the current result. 

• <source>, that contains the name of the analysis, and has attribute type that contains 
the status of the analysis (published or preliminary). 

• <sqrts> contains the collider center-of-mass energy. 

• <CL>: when the likelihood has been extrapolated from a 2D contour, this can be used 
to indicate the CL of the contour thas has been used to extract the covariance matrix 
(usually the 68% or 95% CL contour). 

• <mass> contains the Higgs boson mass considered in the results. If not given, rriu = 
125 GeV is assumed. 
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Note that, if prod="VH" or prod="VVH" is given as attribute to the <expmu> tag or to an 
<ef f > tag, the relative contributions of WH and ZH (for VH) and of WH, ZH and VBF (for 
VVH) will be computed internally assuming an inclusive search, he., for VVH, 

er SM (W) 

eft(V Y) = -- ' ; for A' e (VBF, WH,ZH), (17) 

<Lx=vbf,wh,zh a > 

where the cross sections are evaluated at the Higgs mass given in the <mass> tag using the 
LHC Higgs Cross Section Working Group (HXSWG) results for the 8 TeV LHC [3]. 

An explicit example of well-formed experimental input is 

<expmu dim="2" type="n" decay="ZZ"> 

<experiment>ATLAS</experiments 

<source type="published">HIGG-2013-2l</source> 

<sqrts>7+8</sqrts> 

<mass>125.36</mass> 

<CL>68%</CL> 

<eff axis="x" prod="ggH">l.0</eff> 

<eff axis="y" prod="VVH">l.0</eff> 

<! — (...) “> 

</expmu> 

for the results of the ATLAS H —> ZZ* analysis [36]. The comment <!— (...) —> 

indicates where the likelihood information should be placed. We now present explicitly the 
different possibilities for specifying the likelihood in itself. 

ID interval 

We consider as an example the H —>• bb Tevatron search [52], The following signal 
strength is provided: /r(VH, bb) = 1.59to!72- 

<bestfit>l.59</bestfit> 

<param> 

Uncertainty side="left">-0.72</uncertainty> 

Uncertainty side="right">0.69</uncertainty> 

</param> 

The <bestf it> tag contains the best-fit value for the signal strength. The <uncertainty> 
tag contains the left (negative) and right (positive) lcr errors. If the left and right errors 
are equal in magnitude, the side attribute is not necessary and can be omitted. 



full ID 

We consider as an example the H —>■ bb CMS search [41]. The ID profile likelihood as 
a function of /r(VH, bb) is provided in Ref. [41], see Fig. 2. 

<grid> 

0.00217269 4.30368 

0.00617181 4.26767 


1.99356 3.32783 

1.99756 3.35332 

</grid> 

The digitized likelihood information is stored in the tag <grid>. The first column is 
the signal strength (whose nature is determined by the <eff> tags) while the second 
column is the value of —2 log L. 

2D contour 

We consider as an example the CMS search H —>■ bb [42], The 68 % and 95% CL contours 
are provided in the plane (/i(WH, bb), ^(ZH, bb)). As was explained in Section 2.3, we 
start by fitting the 68 % CL contour assuming that the likelihood follows a bivariate 
normal distribution, and we extract the experimental best-fit point and the inverse of 

the covariance matrix C -1 = 



<bestfit> 

<x>1.123</x> 

<y>0.997</y> 

</bestfit> 

<param> 

<a>l.393</a> 

<b>0.190</b> 

<c>2.217</c> 

</param> 

The tag <bestfit> specifies the location of the best-fit point in the (x,y) plane. The 
tag <param> contains the sub-tags <a>, <b>, and <c>, that parametrize the inverse of 
the covariance matrix in the (x,y) plane. 

full 2D 

We consider as an example the H —» 77 ATLAS search [ 68 ] for which the full likelihood 
information is given in the plane (/i(ggH + ttH, 77, VBF + VH, 77)) [43]. 
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<grid> 

4.26000000e-01 

4.78000000e-01 

5.30000000e-01 


-4.50000000e-01 4.45215260e+01 
-4.50000000e-01 4.19894480e+01 
-4.50000000e-01 3.94115060e+01 


2.87000000e+00 4.45000000e+00 
2.92200000e+00 4.45000000e+00 
2.97400000e+00 4.45000000e+00 


1.84899120e+01 
1.90217220e+01 
1.95605440e+01 


</grid> 


The tag <grid> contains the grid provided by the experimental collaboration. The first 
and second columns are defined by the axis="x" and axis="y" attributes of the <ef f > 
tag, respectively. The third column is the value of —21ogL. 


4.6 User model input 

The user model input, parametrizing the new physics model under consideration, can be 
given either in terms of signal strengths /a(X, Y) directly (defined as in Eq. (6)), or in terms 
of reduced couplings and scale factors (see Section 3). In the latter case, scale factors that 
might be missing in the input are computed, and signal strengths are derived from the scale 
factors. 

The user model input has XML syntax and can be provided as a string or in the form of a 
hie (see the methods readuserinput() and readuserinputfile() in Section 4.2). In this 
section, we present the format that is used in Lilith. 


4.6.1 XML format for signal strengths 

In the signal strengths mode, the basic inputs are the signal strengths defined as in Eq. (6). 
An example of XML input hie for the signal strengths mode is now presented. 


<lilithinput> 

<signalstrengths> 

<mass>125</mass> 


<mu prod="ggH" 
<mu prod="ggH" 
<mu prod="ggH" 
<mu prod="ggH" 


decay="gammagamina">l. 0</mu> 
decay="VV">l.0</mu> 
decay="bb">l.0</mu> 
decay="tautau">l.0</mu> 


<mu prod="VVH" decay="gammagamma">T. 0</mu> 
<mu prod="VVH" decay="VV">l.0</mu> 
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<mu prod="VVH" decay="bb">l.0</mu> 

<mu prod="VVH" decay="tautau">l.0</mu> 

<mu prod="ttH" decay="gammagamina">l. 0</mu> 

<mu prod="ttH" decay="V¥">l.0</mu> 

<mu prod="ttH" decay="bb">l.0</mu> 

<mu prod="ttH" decay="tautau">l.0</mu> 

<!— optionnal: if not given, no decay into invisible —> 
<redxsBR prod="ZH" decay="invisible">0.0</redxsBR> 
<redxsBR prod="VBF" decay="invisible">0.0</redxsBR> 
</signalstrengths> 

</lilithinput> 


• <lilithinput> is the root tag of the XML hie, it defines a Lilith input hie. 

• The <signalstrengths> tag indicates that the user input is given in terms of signal 
strengths. 

• The <mass> tag defines the Higgs boson mass at which the likelihood should be com¬ 
puted. It should be in the [123,128] GeV range. This information is not used in the 
calculations with the current experimental input, where results are only given for a 
hxed Higgs mass. 

• The signal strengths themselves are defined in <mu> tags. Two mandatory arguments 
should be given: 

— The prod attribute can be ggH, WH, ZH, VBF, ttH. For convenience, multi-particle 
attributes have been defined. They are listed in Table 4. 

— The decay attribute can be gammagamma, Zgamma, WW, ZZ, bb, cc, tautau. As 
for the prod attribute, multi-particle labels have been defined and are listed in 
Table 4. 


Note that every <mu> tag can be omitted; in such a case the SM value will be as¬ 
sumed. A warning will furthermore be issued in case of missing Y) for Y G 
( 77 , ZZ*, WW*, bb, tt ) (after resolving multi-particle labels). 


• Finally, there is the possibility to specify an invisible branching ratio in the <redxsBR> 
tag. This is defined as 


redxsBR(A", invisible) 


cjX) 

a SM (X) 


B 


invisible • 


(18) 


As the branching fraction of the SM Higgs boson into invisible particles is very small 
(. B sm (H —> 4 v) = 0.11% at m H = 125 GeV [3]) and cannot be probed at the LHC, 
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Attribute 

shortcut for... 


"VVH" 

"VBF", "WH", "ZH" 

prod 

"VH" 

"WH", "ZH" 


"VV" 

"ZZ", "WW" 

decay 

II II 

"cc", "bb", "tautau" 

"uu" 

"cc" 


"dd" 

"bb", "tautau" 


Table 4: Possible multi-particle attributes for the tag <mu> in the signal strengths mode. 

one usually does not express the results of invisible Higgs searches in terms of signal 
strengths. Invisible decays of the Higgs boson are currently constrained in association 
with two jets from VBF, and in association with a Z boson from ZH production, see 
Table 1. 


Note that the signal strengths for several Higgs states contributing to the signal can be 
defined by specifying an arbitrary number of <signalstrengths> . . </signalstrenths> 
tags in the input. After reading the input, the signal strengths from each individual state 
contributing to the signal will be stored in the attribute user mu, and the sum of the 
signal from the different particules (signal strength per signal strength) will be stored in the 
attribute user_mu_tot (for more details, see Section 4.2). We neglect possible interferences 
between the different states. It can be useful to provide an identifier for each particle. This 
can be achieved with a part attribute to the <signalstrengths> tag. An example of user 
input in terms of signal strengths is stored in userinput/example_mu. xml for the case of 
a single Higgs boson contributing to the signal, and in userinput/examplejnujnultiH.xml 
for the case of two or more particles. 


4.6.2 XML format for reduced couplings 

New physics can be parametrized in terms of scaling factors that can be identified as (or 
derived from) reduced couplings, as was presented in Section 3. In this section we present 
the user input in terms of reduced couplings. Before turning to the format of the user 
input, that also has XML syntax, we comment on the computation of couplings and of signal 
strengths. First of all, as we have seen in Section 3.1, predictions for the Higgs boson can 
be obtained from the reduced couplings Cw , Cz, C t , Cb, C c , and C T appearing in Eq. (9). 
Scaling factors for VBF production and loop-induced processes are function of the C % and 
can be expressed in as Eq. (1 1)—(12) . In the following, we will consider two possible cases: 
that these scaling factors are obtained from leading-order calculations (i.e. tree-level results 
for VBF and one-loop analytical expressions for other processes), or including NLO QCD 
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corrections. The former case will be denoted as LO, the latter one as BEST-QCD. We comment 
on the computations currently implemented in Lilith: 

VBF 

The contribution from the W boson, the one from the Z boson, and the interference 
between them have been obtained from VBFNLO-2 .6.3 [69] for Higgs masses in the 
[123,128] GeV range with (for BEST-QCD mode) and without (for LO mode) NLO QCD 
corrections at the LHC 8 TeV, using the MSTW2008 parton distribution functions [70]. 
The results for cr^y(VBF), crf^f(VBF) and cr^(VBF) as a function of the Higgs 
mass were stored in text files shipped with Lilith and read internally when using 
computereducedcouplingsQ (see Section 4.2). 


ggH 

The contributions from the three heaviest quarks (t, b , c) to the SM cross section are 
taken into account. In the LO mode, we use analytical expressions [58]. In the BEST-QCD 
mode, those have been generated in the [123,128] GeV range with HIGLU [71] at the 
LHC 8 TeV with the MSTW2008 parton distribution functions. 

H -> gg, 77, Z7 

The relevant SM partial widths of these processes (taking into account particles listed in 
Eq. (12)) are obtained from analytical expressions [58] in the LO mode. In the BEST-QCD 
mode, those have been generated in the [123,128] GeV range with HDECAY [72] including 
the available QCD corrections. 

However, the Lagrangian defined in Eq. (9) does not exhaust the possibilities for new physics 
affecting the properties of the Higgs processes. One particularlity interesting case is that 
BSM particles enter the loop-induced processes, such as gg ^ H and H —» 77 . An explicit 
example will be given in Section 6.2. To account for these cases, we allow direct definition 
of scaling factors for the four main loop-induced processes (gg —» H and H —» gg, 77 , Z 7 ), 
i.e. direct definition of C gg H and C 99 i 77 i z 7 . If some or all of the scaling factors are missing 
from the input, they will be computed internally using Eq. (11)—(12), i.e. assuming that only 
SM particles are involved. Finally, note that we use the SM branching ratios provided by 
the LHC HXSWG [3], at the Higgs mass given in the user input, when computing the signal 
strengths (see Eq. ( 8 )). 

The user input file for the Lilith reduced couplings mode has the following structure. 

<lilithinput> 

<reducedcouplings> 

<mass>125</mass> 


<C to="ZZ">1.0</C> 
<C to="WW">1.0</C> 
<C to="tt">l.0</C> 
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<C to="cc">l.0</C> 

<C to="bb">l.0</C> 

<C to="tautau">l.0</C> 

<C to="gammagamma">l.0</C> 

<C to="Zgamma">l.0</C> 

<C to="gg">l,0</C> 

<C to="VBF">l.0</C> 

<precision>BEST-QCD</precision> 

<extraBR> 

<BR to="invisible">0.0</BR> 
<BR to="undetected">0.0</BR> 
</extraBR> 

</reducedcouplings> 

</lilithinput> 


• <lilithinput> is the root tag of the XML hie, it defines a Lilith input hie. 

• The <reducedcouplings> tag is specihc to the reduced couplings mode. This is where 
the reduced couplings are specihed. The correspondence between the XML notation and 
Eq. (9) is given in Table 5. Note the possibility to define common couplings for the 
up-type fermions, down-type fermions, all fermions, and electroweak gauge bosons. 

• The tag <mass> defines the Higgs boson mass at which the likelihood should be com¬ 
puted. The allowed range is [123,128] GeV. This affects the computation of the SM 
branching ratios and partial cross sections and widths as explained above. If it is not 
given, a Higgs mass of 125 GeV is assumed. 

• Regarding the effective coupling to a pair of gluons, NLO corrections affect gluon fusion 
(ggH) and the decay into two gluons (H —> gg ) in a different way. Therefore, scaling 
factors C gg H and C gg can be specihed separately as 

<C to="gg" for="prod">l.0</C> 

<C to="gg" for="decay">l.05</C> 

If for="all" is specihed, the same coupling is assigned to the production and decay 
modes. This is the default behavior if the for attribute is missing. 

• CP violation was presented in Section 3.2. In the LO mode, the fermionic couplings 
C t , Cb , C c , C T can be given a real and an imaginary component. For the top quark for 
instance, this can be specihed as 
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to attribute 

"tt" 

"bb" 

"cc" 

"tautau" 

"WW" 

"ZZ" 

corresponds to 

c t 

c b 

C c 

C T 

Cw 



to attribute 

"gammagamma" 

"Zgamma" 

"f 

for="prod" 

5g" 

for="decay" 

"VBF" 

corresponds to 

c 

Cz'y 

CggH 

c 99 

Cvbf 


to attribute 

shortcut for... 

> 

> 

"WW", "ZZ" 

II -p -p 1! 

"tt", "cc", "bb", "tautau" 

"uu" 

"tt","cc" 

"dd" 

"bb", "tautau" 


Table 5: Upper and middle tables: correspondence between the XML notation and that of 
Eq. (9). Bottom table: common reduced couplings definitions. 

<C to="tt" part="re">0.8</C> 

<C to="tt" part="im">0.2</C> 

If part="re|im" is not specified, the coupling is assumed to be purely real. In the 
BEST-QCD mode, only the real part of the coupling is taken into account. 

• The <precision> tag contains either BEST-QCD or LO. If not specified, or wrongly 
spelled, the BEST-QCD mode is the default mode. 

• The <extraBR> tag contains the declaration of the invisible or undetected branching 
ratios (see Section 3). 


As in the case of input in terms of signal strengths, several tags <reducedcouplings> . . 
</reducedcouplings> can be defined, corresponding to the case where several Higgs states 
contribute to the observed signal around 125 GeV. These particles can also be given a name 
with a part attribute to the <reducedcouplings> tag. An example of user input in terms 
of reduced couplings is stored in userinput/example_couplings .xml for the case of a single 
Higgs boson contributing to the signal, and in userinput/example_couplings_multiH. xml 
for the case of two or more Higgs states. 
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4.7 Output 


When using the API, all relevant information can be accessed from the public attributes of the 
class Lilith presented in Section 4.2. The user can manipulate and store this information in 
any way he/she wants. However, having standardized output formats is important in order to 
interface Lilith with other programs. We provide three output methods: writecouplings(), 
writesignalstrengths(), and writeresults() (for details on how to call these methods, see 
Section 4.2). The first two methods write the reduced coupling information (if available) and 
the signal strength information, respectively, in a hie. The formats specified in Section 4.6.1 
and 4.6.2 above are used, such that these output hies can also be used as input to a subsequent 
call to Lilith. In the command-line interface, these two methods can be called with the 
options -c or —couplings, and -m or —mu, respectively (for more details, see Section 4.3). 

The last method, writeresults(), is used to store the results after evaluation of the likelihood. 
Depending on the second argument, the output will be written in XML format (if slha=False ) or 
in SLHA-like format (if slha=True). By default the output hie is in XML format for consistency 
with what is used otherwise in Lilith; an SLHA-like output was added as it is widely used 
in BSM phenomenology. This method can be called with the option -r or —results in the 
command-line interface. 

We start with the description of the XML format. Its root tag is <lilithresults>. The results 
from each analysis are given within an <analysis> tag having two attributes: experiment 
and source, corresponding to the value of the tags <experiment> and <source> read from 
the experimental input, if present; otherwise the attribute is given an empty value Each 
<analysis> tag contains a <1> tag, whose value is —21ogL from this experimental result 
alone, and an <expmu> tag that follows the syntax of the experimental input (see Section 4.5) 
except that it only provides information on the efficiencies (in <eff> tags). 

Finally, <lilithresults> contains tags summarizing the information: <ltot>, whose value 
is the sum of the —21ogL values from all experimental results, and <expmdf>, the number 
of measurements. We also provide the version of Lilith in a <lilithversion> tag and the 
version of the experimental database of results in a <dbversion> tag. To summarize, we give 
a complete example of output considering only two experimental results. 

<lilithresults> 

<lilithversion>l.l</lilithversion> 

<dbversion>15.02</dbversion> 

<analysis experiment="ATLAS" source="HIGG-2013-08"> 

<expmu decay="gammagamma" dim="2" type="n"> 

<eff axis="x" prod="ggH">l.0</eff> 

<eff axis="y" prod="VBF">l.0</eff> 

</expmu> 

<1>0.9</l> 
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</analysis> 


<analysis experiment="CMS" source="CMS-HIG-13-030"> 

<expmu decay="invisible" dim="l" type="f"> 

<eff prod="VBF">0.79</eff> 

<eff prod="ZH">0.21</eff> 

</expmu> 

< 1 > 1 . 1 </ 1 > 

</analysis> 

<ltot>2.0</ltot> 

<exp_ndf>3</exp_ndf> 

</lilithresults> 

The output in SLHA-like format is much more basic. Currently, we only provide the total 
—21ogL value as well as the number of measurements. Taking the same example as above, 
the output would read: 

BLOCK LilithResults 

0 2.0 # -2*LogL 

13 # exp_ndf 

The SLHA structure, where each element of a block is identified by a set of numbers, makes 
it more complicated to integrate all necessary information in a well-structured way. For this 
reason we recommend the XML format. Extensions of the SLHA format will be considered in 
the future depending on the needs of the users of Lilith. 


5 Validation 

Having explained how to use Lilith in the previous section, we now turn to the validation 
of the likelihood derived from the experimental input shipped with the code. We begin 
by discussing the validity of the bivariate normal distribution as an approximation to the 
2D likelihood functions in the signal strength planes (^(X, Y'), /i(X',Y')). The use of this 
approximation is necessary whenever only contours of constant likelihood are provided instead 
of the full information. Several coupling fits from the ATLAS and CMS collaboration are 
then reproduced. The results from Lilith are compared to the official ones to assess the 
validity of the likelihood used in Lilith. 
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Figure 3: Reconstruction of the experimental likelihood from a bivariate normal approxi¬ 
mation for the ATLAS WW* search [47] (left) CMS 77 search [34] (right). The filled dark 
and light grey contours show the 68 % and 95% CL experimental contours while the red and 
orange solid lines show the reconstructed likelihood contours. The blue diamond and the 
black star indicates the experimental and reconstructed best-fit points, respectively. 


5.1 Reconstruction of the experimental likelihoods 

In a signal strength plane (n(X,Y), n(X 'an approximation to the likelihood function 
can be obtained assuming that the measurements follow a bivariate normal distribution, as 
explained in Section 2.3. Using the 68 % CL contour provided by the experimental collabo¬ 
ration, we reconstruct the shape of the likelihood and compare the location of the best-fit 
point as well as the 68 % and 95% CL contours with what is provided by ATLAS or CMS. 

Two examples are shown in Fig. 3: the reconstruction of the likelihood for the ATLAS 
WW* [47] and the CMS 77 [34] final states.' In both cases we observe an excellent agreement 
between the reconstructed likelihood and the official result. The 68 % CL regions are perfectly 
reproduced and the reconstructed best-fit points are very close to the experimental ones. 
The extrapolation towards the 95% CL regions also shows very good agreement. We find 
equally good agreements with all other decay modes (with the exception of H —y ZZ*), 
and we conclude that the Gaussian distribution is a very good approximation to the true 
distribution. 

The largest deviations from the normal approximation are expected to occur for final states 
with low statistics since the counting of the events, that follows the Poisson distribution, has 
not yet entered the Gaussian regime. In particular, this is the case for the ZZ* channel. In 

'Additional validation materials for the reconstruction of the likelihood can be found at [67]. 
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Figure 4: Reconstruction of the experimental likelihood from a bivariate normal approxima¬ 
tion for the ATLAS [36] (left) and CMS [53] (right) H —>• ZZ* searches. The filled dark 
and light grey contours show the 68% and 95% CL experimental contours while the red and 
orange solid lines show the reconstructed likelihood contours. The blue diamond and the 
black star indicates the experimental and reconstructed best-fit points, respectively. 


Fig. 4, we show the comparison between the Lilith reconstructed likelihood in the ZZ* final 
state and the corresponding ATLAS and CMS ones. 

As can be seen, the deviation of the ATLAS likelihood from the bivariate normal approxi¬ 
mation can be substantial. In the positive region of the plane (the one that is relevant), the 
approximation holds well near the best-fit point. However, going away from it the recon¬ 
structed shape fails to reproduce the ATLAS 95% CL contour at large /i(VBF + VH, ZZ*). 
Due to non-Gaussian effects, the reconstructed best-fit point is quite distant to the exper¬ 
imental one. For the CMS case, the approximation holds to a better approximation. The 
reconstructed best-fit point is very close to the experimental one and the shape of the 95% CL 
contour is very well reproduced although a small shift in the /i(ggH + ttH, ZZ*) direction is 
observed. As will be argued in Section 7, provision of the full likelihood information would 
yield a significant improvement over the normal approximation in such cases. 


5.2 Comparison to Higgs coupling fits from ATLAS and CMS 

In order to validate the approximate Higgs likelihood used in Lilith, we attempt to repro¬ 
duce coupling fit results from combination notes from ATLAS [4] and CMS [5]. Note that 
while the CMS combination [5] makes use of the final Run I results, a number of analyses 
considered in the combination of the ATLAS results given in Ref. [4] have been updated 
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Figure 5: (C 7 , C g ) (left) and (Cy, Cf) (right) fits using data from the ATLAS combination [4], 
The red and orange filled surfaces correspond to the 68% and 95% CL regions obtained by 
the ATLAS collaboration while the corresponding dashed lines show the Lilith results. The 
black star indicates the position of the Lilith best-fit point, the blue diamond is the ATLAS 
best-fit point and the white circle shows the SM prediction. 


since then. The final, legacy combination of the Higgs measurements from ATLAS at Run 1 
has not yet been released. For this reason, the lists of recommended experimental results 
data/CMS-HIG-14-009. list and data/latestCMS. list are identical, as of February 2015, 
while data/ATLAS-C0NF-2014-009. list and latestATLAS. list differ. 

First, results for two benchmark scenarios proposed by the LHC HXSWG in [73] are pre¬ 
sented. In the first scenario, SM-like tree-level couplings are assumed (i.e., all C t — 1 
in Eq. (9)) but two scaling factors are introduced: C 7 = C 77 (scaling H —> 77 ), and 
C g = CggH = C gg (scaling ggH production and H —> gg ). In the second benchmark sce¬ 
nario, two reduced couplings are introduced: Cy = Cw = Cz, for the coupling of the Higgs 
boson to a pair massive vector bosons, and Cf = C t = Cb = C c = C T , a universal coupling to 
fermions. In this case, the effective coupling to gluons is simply Cp, while C 77 is a function of 
both Cy and Cp that was obtained taking into account QCD corrections (see Section 4.6.2). 

Let us first discuss the results from ATLAS, obtained using the list of experimental results 
data/ATLAS-C0NF-2014-009. list. Results of the two fits are presented in Fig. 5. In both 
scenarios, very good agreement is observed between the results from ATLAS and the ones 
obtained with Lilith. Both the reconstructed best-fit point and contours reproduce very 
well the ATLAS results. The most significant deviation is a slight deformation of the 95% CL 
region in the (C' 7 , C g ) plane. The corresponding results for CMS are shown in Fig. 6 . CMS 
results are well reproduced with Lilith, even for the contour at 99.7% CL. Slight shifts of the 
best-fit points and minor deformations of the contours are observed. The overall agreement 
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Figure 6: (C^,C g ) (left) and ( Cv,Cf ) (right) fits using data from the CMS combination [5]. 
The red, orange and yellow filled surfaces correspond to the 68%, 95% and 99.7% CL regions 
obtained by the CMS collaboration while the corresponding dashed lines show the Lilith 
results. The black star indicates the position of the Lilith best-fit point, the blue diamond 
is the CMS best-fit point and the white circle shows the SM prediction. 


is nevertheless very good. 

Let us move on to the 3-parameter fit (Cw, Cz, Cp). As in the (Cy, Cp) benchmark scenario 
discussed above a universal coupling to fermions is introduced, but instead of a single coupling 
to vector boson one defines separately the reduced coupling to W bosons, Cw, and to Z 
bosons, Cz Defining Cwz = Cw/Cz, a direct test of custodial symmetry can be performed 
using the Higgs measurements alone. The 1-dimensional likelihood profile for Cwz is shown 
in Fig. 7 for both the ATLAS and CMS combination. 

Although the ATLAS result is almost perfectly reproduced, a significant discrepancy is ob¬ 
served in the case of CMS for Cwz > 1- This does not come as a surprise: several experimen¬ 
tal results were considered in the (/i(ggH + ttH, Y), /r(VBF + VH,F)) plane. The breaking 
of VBF + VH into the individual production modes VBF, WH and ZH (assumed to be inclu¬ 
sive, see Eq. (17)) becomes relevant for C w % Cz- Moreover, ATLAS results make use of the 
full numerical likelihood grids that were provided in Refs. [43-45] while the bivariate normal 
approximation is used in the case of CMS. Thus, constraints on models in which Cw % Cz 
should be interpreted with care given the experimental information being used as input to 
Lilith. 

Finally, we present the result of a 3-parameter fit (C 7 , C g . £>i nv i S ibi e ) in terms of the ID profile 
likelihood of ^invisible in Fig. 8. A very good agreement is observed in ATLAS, and in CMS 
for moderate values of i3 invisi bi e - As explained in Section 3.1, the presence of a branching 
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Figure 7: ID likelihood profiles of C-wz = C-w/Cz from a ( C Wl C Zl Cp) fit to the ATLAS [4] 
(left) and CMS [5] (right) data and comparison to the official results. The ATLAS fit considers 
both signs for the Higgs-fermion-fermion coupling and furthermore defines Cfz = Cf/Cz- 
The results are given for both signs of Cfz- 



Figure 8: ID likelihood profiles of Lwisibie from a (C 7 , C g . invisible) fit an d comparison to the 
ATLAS [4] (left) and CMS [5] (right) results. 
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Figure 9: 1-dimensional likelihood profiles of 
benchmark scenario (CV,Cp). 



V (left) and Cf (right) from a global fit of the 


ratio into invisible particles is constrained by direct searches for invisible decays of the Higgs 
boson, and also by every Higgs search since it modifies the total Higgs width and therefore 
scales all signal strengths collectively. 


6 Examples of applications 

Having validated the Higgs likelihood in Lilith from results obtained by the ATLAS and 
CMS collaborations, we now turn to deriving constraints on specific new physics scenarios 
using the latest LHC results (as of February 2015) present in data/latest. list. The Python 
routines used to obtain these results are available in the folder examples/python and will be 
described shortly. 

6.1 Reduced coupling determination 

As a first illustration of the use of Lilith, constraints on the benchmark scenario ( Cy,Cp ) 
introduced in Section 5.2 are derived. The right panels of Figs. 5 and 6 show results on this 
scenario in the 2D plane (CV, Cf) using only ATLAS or CMS results. Here we combine the 
latest ATLAS and CMS results and derive ID profile likelihood constraints on Cy and Cf- 
Results are shown in Fig. 9. 

The Python routine used to obtain this result is CVCF_ldprof ile .py. ft can be executed 
from the Lilith-1.1/ folder with the command line 
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i [python examples/python/CVCF_ldprofile.py 

This example uses the class Minuit of the library iminuit [74], a Python implementation of 
the MINUIT [75] minimization library, in order to minimize —21ogL and derive the ID profile 
around the minimum. Moreover, matplotlib [76] is used to produce the resulting figures. 
Below, we describe parts of the routine. 

After having instantiated the Lilith class and read the experimental data with 

1 lilithcalc = lilith.Lilithfverbose, timer) 

2 lilithcalc.readexpinput(myexpinput) 

a function getL returning —21ogL for a given (Cy, Cp) point is defined 

3 def getL(CV, CF): 

4 myXML_user_input = usrXMLinput(mh=mh, CV=CV, CF=CF, precision=precision) 

5 lilithcalc.computelikelihood(userinput=myXML_user_input) 

e return lilithcalc .1 

where the function usrXMLinput creates a XML user input string from Cy and Cp, for a given 
precision mode precision. 

An object m of the class Minuit is then created 

7 |m = MinuitfgetL, CV=1, limit_CV=(0,3), CF=1, limit_CF=(0,3)) 

where the initial point of the minimization and the range of parameters are defined. The 
function m.mnprofile is then called 

s |xV, yV, rV = m.mnprofile("CV", bins=300, bound=(0., 2), subtract_min=True) 

and returns the ID likelihood profile A(—2 log L(Cy)) = —2log(L(Cy)/L(best fit)) for a 
given range and number of points. 8 Without the option substract_min=True, the “absolute” 
likelihood —2 log L(CV) would be returned instead. 

The parameter range in which A(—2 log L(Cy)) < 1 (4) defines the 68% (95%) CL intervals 
of Cy. The constraints on Cp are derived in the same way, and all results are plotted and 
stored in results/CVCF_ldprof ile. pdf. In this scenario, the best-fit point is obtained for 
Cy = 1.05 and Cf = 1.02. In ID, i.e. profiling over the other parameter, the 68% (95%) CL 
intervals read C v = [1.00,1.09] ([0.96,1.13]) and C F = [0.92,1.12] ([0.83,1.22]). 

We also provide example on how to derive constraints and produce figures for a 2D parameter 
space. The left panel of Fig. 10 presents the 2D contraints obtained from a global fit of the 
(Cj,C g ) model presented above. The corresponding Python routine is CgammaCg_ 2 d.py. It 
can be executed from the Lilith-1.1/ folder with the command line 

8 In general for a function —21og L({Ci},{Cj}), where {Ci},{Cj} can be sets of parameters, the profile 
likelihood — 21ogL({Ci}) is obtained by minimizing the full function with respect to {Cj} for a given {Cj}- 
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Figure 10: Contraints on (C 7 , C g ) (left) and (CV,Cf) (right) from a global fit to the Higgs 
data. The red, orange and yellow filled surfaces correspond to the 68%, 95% and 99.7% CL 
regions. The black star shows the position of the best-fit point. 


i [python examples/python/CgammaCg_2d.py 

A scan of the (C 7 , C g ) parameter space is performed, and results are stored in the hie 
results/CgammaCg_2d.out in the form 


0.04040 0.00000 
0.04040 0.02020 
0.04040 0.04040 
0.04040 0.06061 


119.05462 

119.00658 

118.86261 

118.62314 


where the first, second and third columns contain the values of C 7 , C g and —21ogL(C' 7 , C g ), 
respectively. The 68%, 95%, 99.7% CL regions in the (Cj,C g ) plane then corresponds 
to A(—21og L(C^,C g )) < 2.3,5.99,11.83, respectively. This identification is performed by 

matplotlib with 

i | ax .contour (xi,yi,Z, [2.3,5.99,11.83]) 

where xi, yi and Z are list of points defining the grid in the C 7 and C g directions, and 
the corresponding A(—2 logL(C' 7 , C g )) value, respectively. The results are displayed and 
stored in results/CganunaCg_2d.pdf. For completeness, the 2D constraints on the (CV,C(f) 
benchmark scenario, using the latest LHC measurements, are also presented in the right 
panel of Fig. 10. They have been derived in the same way. 
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6.2 Higgs constraints on superpartners of the tau lepton 

Supersymmetric scalar partners of the tau leptons, known as staus, can have substantial 
contribution to the H —y 77 decay rate if they are light and have a large mixing [77,78]. 
Constraints on the parameters controlling this new contribution can therefore be obtained 
from the Higgs precision measurements. Here, we consider the Minimal Supersymmetric 
Standard Model (MSSM) and assume that the only deviation from a SM-like Higgs behavior 
comes from the contribution of staus to the loop-induced process H —* 77 . More precisely, 
it is assumed that the supersymmetric partners of the Higgs boson and of the remaining 
fermions are decoupled, that the second Higgs doublet is phenomenologically irrelevant and 
that a Higgs mass of 125 GeV can be obtained for any point of the analysis. In this case, 
the contribution from staus to the H —y 77 decay width is parametrized by the two phys¬ 
ical masses rrif 1 and rrif 2 (with rrif 1 < raf 2 ), the mixing angle Of and the ratio of vacuum 
expectation values for the two Higgs doublets, tan j3. The corresponding amplitude at LO 
reads [58,79] 




= E 

i=l,2 


9hi 


j, Of, tan /3) 


m% 


A h 




(19) 


where the sum runs over the two stau mass-eigenstates, Aq is a loop form factor and g h f z f 
is the Higgs stau stau coupling. 

The effective Higgs- 7-7 reduced coupling can therefore be expressed as 


C-y {rrif 1 , m,f 2 , Of, tan f3) 


•MffTY + ^if 77 ( m fi) 

m,f 2 , Of, tan j3) 


A/fSM 
JVI H 77 



( 20 ) 


Note that the SM amplitude appears both in the numerator and denominator of 

Eq. (20) since SM tree-level couplings are assumed. 

Fixing tan/3 = 10 and the mass of the lightest stau rrif 1 = 85,100 GeV, we show constraints 
in the plane (rrif 2 , Of) in Fig. 11. For Of = 7 t/4 , the 2-dimensional 95% CL upper limit on 
rrif 2 reads rrif 2 < 360 (460) GeV for m fl = 85 (100) GeV. More generally, the upper limit on 
rrif 2 becomes weaker as m fl is increased. 

The corresponding Python code is stau_gammagamina. py. It can be executed by typing the 
following command line to the shell from the Lilith-1.1/ folder: 

1 I python examples/python/stau_gammagamma.py 


The routine works as follows. Functions returning C 1 according to Eq. (20) and —2 log L(C 7 ) 
are defined. Since tan / 3 and rrif 1 are fixed, a 2 -dimensional grid scan is then performed over 
the two remaining parameters: for each couple (■ rrif 2 ,0f ), the corresponding A(—21ogL) is 
obtained. The 2-dimensional 68 %, 95%, 99.7% CL regions in the plane [rrif 2 ,0f) are obtained 
with A(—21ogL) < 2.3,5.99,11.83, respectively. 
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Figure 11: Constraints on the staus masses and mixing angle in the ( m~ 2 ,0 ~) plane for 
m~ = 85 GeV (left) and m~ = 100 GeV (right) and tan /3 = 10. The red, orange and yellow 
filled surfaces correspond to the allowed 68%, 95% and 99.7% CL regions, respectively. 


Note that direct searches from LEP [80] and vacuum metastability condition [81] impose 
further constraints on this scenario. Moreover, this simplified SUSY scenario could easily be 
generalized, e.g. by taking into account H —» XiXv Light staus are especially relevant in 
the case where is light in order to have a viable neutralino dark matter candidate in the 
MSSM (see, e.g., Ref. [82] and references therein). 


7 Prospects for Run II of the LHC 

As we discussed in Section 2, approximations necessarily need to be made when combining 
signal strength results from several categories or several searches, making it necessary to 
validate the approach. In Section 5, we have shown that we reproduce well the results of 
coupling fits from ATLAS and CMS (separately). However, it is clear that the situation 
will change as more statistics will be collected at Run II of the LHC. Indeed, systematic 
uncertainties will then dominate over statistical uncertainties in the majority of the channels. 
Missing correlations between systematic uncertainties (both theoretical and experimental) 
will thus become a more pressing issue. Moreover, more combinations for production and 
decay of the Higgs boson, (X,Y), will be determined with a good precision. This will spoil 
the simple interpretation we have for a number of the results we currently use, in particular 
for results given in the plane (/i(ggH + ttH, Y), ^(VBF + VH, Y)). In this section we recall 
the main limitations when using the information currently provided by the ATLAS and CMS 
collaborations for constructing a likelihood. We also discuss new ways of presenting the LHC 
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Higgs results in order to be able to construct a good approximation to the Higgs likelihood at 
Run II of the LHC. This section is partly based on the note “On the presentation of the LHC 
Higgs Results” [83] that was put forth by a collaboration of theorists and experimentalists 
with the aim to maximize the impact of the LHC Higgs results and their utility to the whole 
high-energy physics community. 

First of all, in most cases only contours of constant likelihood (at least the 68 % CL interval 
or contour, sometimes contours at 95% CL) are provided by ATLAS and CMS. This makes 
it necessary to extrapolate the likelihood assuming, most naturally, a Gaussian shape. When 
using a given contour to extrapolate the likelihood, the validity of this approximation can 
be tested from a comparison of the position of the best-fit point and from contours provided 
by the experimental collaboration. This was done in Section 5.1, where we concluded that 
the reconstruction is generally very good, although in some cases asymmetrical effects are 
washed out (see, e.g., Fig. 4). However, in all cases it induces an unnecessary source of error. 
The ATLAS and CMS collaborations initiated some efforts during Run I to provide the full 
likelihood information in ID and 2D planes (see the right panel of Fig. 2 and Refs. [43-45]). 
We strongly hope that this will become standard practice during Run II of the LHC, and 
that the information will systematically be provided in numerical form. 

Another issue is the dependence of the results on the assumed Higgs boson mass mu- Cur¬ 
rently, we use results given at a fixed Higgs mass. As not all results are provided at the 
same m#, a slight inconsistency is introduced in the combination of the different results 
(the assumed Higgs mass varies within a few hundreds of MeV). Official combination notes 
allow us to get rid of this inconsistency, as all results are therein given at the same Higgs 
mass. However, the dependence of the experimental results on the Higgs mass can be very 
important for the high-resolution channels, that target decay of the Higgs boson into charged 
leptons and photons (such as H —* ZZ* —* 4£ and H —> 77 ). Thus, it would be highly 
desirable to have access to mass-dependent likelihood results. 

Current results are presented in 1 - or 2 -dimensional projections, often corresponding to the 
combination of production modes (in 2D, typically ggH + ttH and VBF + VH). As we dis¬ 
cussed above, this becomes a limitation as measurements get more precise, in which case 
we would like to investigate deviations in all of the five production modes separately. For 
such reasons, a total breakdown of the signal strength measurements in terms of the five 
Higgs production modes (ggH, VBF, WH, ZH, ttH) would be a considerable step forward 
regarding the interpretation of the LHC Higgs results. We would therefore like to advocate 
the experimental collaborations to provide the likelihood as a function of the Higgs mass and 
a full set of production modes, that is to say, in the 

(m H , MggH, n MttH, n MV BF, V), MZH, V), M WH, V)) ( 21 ) 

parameter space for each final state Y. For some final states, all five production modes are 
certainly not constrained by the experimental searches and only lower dimensional projections 
of this space would be relevant. 

This would solve most of the limitations currently faced, with the notable exception of cor- 
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relations between the measurements of different decay modes. For instance, theoretical un¬ 
certainties on gluon fusion production affect both the 77 and ZZ * final states. Recently, an 
interesting proposal was made in this direction in Ref. [84]. Provided experimental collabo¬ 
rations publish likelihoods that are not profiled over a set of theoretical nuisance parameters 
of interest, but instead given for a fixed scenario, it is possible to build a “recoupled” like¬ 
lihood incorporating these uncertainties at the later stage. This has the advantage of not 
being restricted to the Gaussian approximation. It would certainly be of great interest if 
the information in the 2D plane (/x(ggH + ttH, Y), /i(VBF + VH,1')), or even better in the 
possibly 6 D plane discussed above, could be given without profiling over the theoretical un¬ 
certainties on the Higgs signal. With the method presented in Ref. [84], one could then fully 
correlate the theoretical uncertainties between the different channels and experiments, and 
modify these uncertainties compared to what is done in ATLAS and CMS if desired. 


8 Conclusions 

Crucial information on the origin of the electroweak symmetry breaking—and, more generally, 
on the presence of BSM physics around the electroweak scale—can be obtained from the 
study of the properties of the Higgs boson with mass around 125 GeV at the LHC. The 
presentation of the experimental results in terms of signal strengths makes it possible to 
combine all measurements and perform a global fit to the properties of the Higgs boson. The 
results of such fits can be used to discriminate between models where the structure of the 
couplings to SM particles is as in the SM. 

However, using all available information from the ATLAS and CMS collaborations to con¬ 
struct a likelihood is a non-trivial task. Indeed, there is a wealth of experimental searches, 
from which the necessary information often needs to be extracted and put in numerical form. 
Care should also be taken in order to include all available correlations between systematic 
uncertainties. To this aim, we provide a new public tool, Lilith. Lilith is a library written 
in Python, and for which we provide an API as well as a command-line interface and a basic 
interface to C and C++/R00T. The experimental results are read from a database in XML format 
that is shipped with the code and which is easy to modify and extend. Lilith uses as a 
primary input results in which the fundamental production and decay modes are unfolded 
from experimental categories. 

New physics can be parametrized in terms of reduced couplings, or signal strengths directly, 
which are given as input to Lilith in XML format. If needed, scaling factors for the loop- 
induced processes and VBF production are computed taking into account QCD corrections. 
CP-violating Higgs couplings can also be given as input to Lilith. The likelihood is evaluated 
from a set of experimental results and given as output; detailed results can moreover be stored 
in XML or SLHA-likc format. For convenience, Lilith is provided with several applications of 
the code where constraints on effective or explicit models of new physics are derived. They 
include extensive comments. 
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The Higgs likelihood of Lilith obtained from the latest measurements at the LHC has 
been thoroughly validated against ATLAS and CMS results and can be used to constrain 
new physics. Future measurements at Run II of the LHC will, however, call for new ways 
of presenting results in order to derive a good approximation to the Higgs likelihood. In 
particular, further disentanglement of the different production and decay modes will become 
necessary. Moreover, correlations between systematic uncertainties, and in particular the 
treatment of theoretical uncertainties, will become a more pressing issue. The structure 
of the code is such that Lilith can easily be adapted to handle extended signal strength 
information. 
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